###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##
###*### ESTIMATE LOCAL-SPECIFIC PARAMETERS                #*###*###*###*###*###*###*###*##
###*### USING PREVIOUSLY ESTIMATED NAC PARAMETERS         #*###*###*###*###*###*###*###*##
###*### LOCATION 4: HIGH I - HIGH P                       #*###*###*###*###*###*###*###*##
###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*###*##

#' FROM THE PAPER:
#' Child Care Markets, Parental Labor Supply, and Child Development
#' BY:
#' Samuel Berlinski, Maria Marta Ferreyra, Luca Flabbi and Juan David Martin (2022)

# ___________________________________________________________________________________ ####
# START ####

library(lubridate)
library(snowfall)
library(readxl)

rm(list = ls())
gc()

## Directory
folderpath <- '~/EstimationCodeLocation/'
setwd(folderpath)

# ___________________________________________________________________________________ ####
# INPUTS ####

## Functions
source('functions.R')
source('simul_market.R')

## Sim data
load('SimData_L4.RData')

## Moments
moments     <- as.data.frame(read_excel('Location-SpecificMoments.xlsx', sheet = 'L4'))
sam_moments <- round(moments$Value, 2)
weights     <- 1/moments$SE#moments$Weight

# ___________________________________________________________________________________ ####
# SETUP ####

## Market given settings
J     <- nrow(data)
NL    <- 5
ratio <- 1.0593
minw  <- 5.15*ratio
M     <- 40*NL/J/(0.3 * 1.1117080)
Fixed <- c(128.90, 322.25)*ratio
Cn    <- minw
Cr    <- 0
I     <- minw*8
subs  <- 0.35
Nrel  <- round(J*0.40*0.9196207)

## Observed heterogeneity
q0  <- data$q0
ed  <- data$edu
edf <- data$edf
wf  <- data$wf
sm  <- data$single
To  <- 16 + 2*sm

## Dad's wage implied random term
ef  <- log(wf)      
mwf <- mean(ef[sm == 0])
swf <- sd(ef[sm == 0])
ef  <- (ef - mwf)/swf
ef[wf == 0] <- 0

## Random draws
set.seed(123)
ev   <- matrix(rnorm(3*J), J, 3)
ew   <- matrix(rnorm(J),   J, 1)
ec   <- rowMeans(apply(matrix(rnorm(NL*10000), NL, 10000), 2, sort))
ec   <- (ec - mean(ec))/sd(ec)
reli <- sample(1:J, size = Nrel, replace = F)

## Standarization for wages draws
ew[ed == 1] <- (ew[ed == 1] - mean(ew[ed == 1]))/sd(ew[ed == 1])
ew[ed == 2] <- (ew[ed == 2] - mean(ew[ed == 2]))/sd(ew[ed == 2])
ew[ed == 3] <- (ew[ed == 3] - mean(ew[ed == 3]))/sd(ew[ed == 3])

# ___________________________________________________________________________________ ####
# ESTIMATE LOCAL PARAMETERS ####

## Nac parameters & local parameters
load('Estimates_01.RData')
pp1 <- c(pp1, 'dwm' = ratio, 'dws' = ratio, 'dcm' = ratio, 'dcs' = ratio, 'drho' = 1)

## Estimate by minimizing the simulated moments function
keep <- c(objects(), "keep")
sfInit(parallel = TRUE, cpus = 3)
sfExport(list = keep, local = FALSE)
n   <- c('dwm', 'dws', 'drho', 'dcm', 'dcs')
ans <- optim(pp1[n], fn = crit_funn, parameters = pp1, 
             nampar  = n, sample_moments = sam_moments, 
             weights = weights, method = 'Nelder-Mead', 
             control = list(reltol = 1e-5, maxit = 1000))
sfStop()

## Save the new vector of parameter estimates as the optimal solution to eq. (26)
pp1[n] <- ans$par
pp1    <- conv(pp1)
pp1    <- constraints(pp1)
pp1    <- conv.inv(pp1, local = T)
save(pp1, file = 'EstimatesL4_01.RData')

# ___________________________________________________________________________________ ####
# END ####